Method and system for imaging the dynamics of scattering medium

ABSTRACT

A method and system for imaging the dynamics of a scattering medium ( 116 ) is provided. The method and system generates contrast and resolution enhanced images of dynamic properties of a medium having a temporal signature by using time series analysis methods on a time series of collected data or time series of images to extract and isolate dynamic properties of the medium ( 116 ).

This application claims the benefit under 35 U.S.C. §120 of prior U.S.Provisional Patent Application Ser. No. 60/153,926 filed Sep. 14, 1999,entitled DYNAMIC TOMOGRAPHY IN A SCATTERING MEDIUM and 60/154,099 filedSep. 15, 1999, entitled DYNAMIC TOMOGRAPHY IN A SCATTERING MEDIUM.

This application is related to copending application Ser. No.PCT/US00/25155, filed on the same date as this application, entitled“SYSTEM AND METHOD FOR TOMOGRAPHIC IMAGING OF DYNAMIC PROPERTIES OF ASCATTERING MEDIUM” by inventors R. Barbour and C. Schmitz and is herebyincorporated by reference (hereinafter the “Barbour 4147PC1application”). The counterpart U.S. patent application is applicationSer. No. 10/088,254 filed Mar. 14, 2002.

This application is also related to copending application Ser. No.PCT/US00/25156, filed on the same date as this application, entitled“IMAGING OF SCATTERING MEDIA USING RELATIVE DETECTOR VALUES” by inventorR. Barbour and is hereby incorporated by reference (hereinafter the“Barbour 4149PC2 application”). The counterpart U.S. patent applicationis application Ser. No. 10/088,192, filed Mar. 14, 2002.

This invention was made with U.S. Government support under contractnumber CA-RO166184-02A, awarded by the National Cancer Institute. TheU.S. Government has certain rights in the invention.

FIELD OF THE INVENTION

The invention relates to the field of imaging in a scattering medium,and in particular to methods of imaging the dynamics of a scatteringmedium.

BACKGROUND

Imaging in a scattering medium relates generally to the methods andtechniques for generating an image of the internal properties of ascattering medium based on the detection of scattered energy.

Many systems and techniques have been developed for imaging ofscattering media. A typical system for imaging based on scattered energydetection includes a source for directing energy into a target mediumand a plurality of detectors for measuring the scattered energy exitingthe target medium at various locations with respect to the source. Basedon the measured energy exiting the target medium, it is possible toreconstruct an image representation of the cross-sectional scattering,absorption or other properties of the target medium. The values of theabsorption and scattering properties of the medium can vary depending onthe wavelength and types of energy employed as an imaging source. Thesevalues are also frequently spatially varying. These techniques permitthe use of types of energy and wavelengths, such as near infrared lightenergy, that are not suitable for projection imaging techniques, such asx-ray imaging. Thus these techniques have great potential for detectingand imaging properties of media, such as human tissue, that can not berevealed using energy sources commonly employed in projection imagingmethods.

Exemplary methods and systems for imaging of scattering media aredisclosed in Barbour et al., U.S. Pat. No. 5,137,355, entitled “Methodof Imaging a Random Medium,” (hereinafter the “Barbour '355 patent”),Barbour, U.S. Pat. No. 6,081,322, entitled “NIR Clinical Opti-ScanSystem,” (hereinafter the “Barbour '322 patent”), the Barbour 4147PC1application, the Barbour 4149PC1 application and the Barbour 4147PC2application.

As can readily be appreciated, there are many instances where use ofthese techniques are highly desirable. For example, one flourishingapplication area is in the field of optical tomography. Opticaltomography typically uses near infrared (NIR) energy as an imagingsource. Contrary to imaging methods relying on the use of ionizingradiation and/or toxic/radioactive contrast agents, NIR opticaltomographic imaging methods bear no risk of causing harm to the patient.The dose of optical intensity used remains far below the threshold ofthermal damage and is therefore safe. In the regime ofwavelength/intensity/power used, there are no effects on the tissue thataccumulate with increasing light energy dose due to over-all irradiationtime.

Other favorable attributes of optical tomography include the use oflow-cost, potentially portable devices that employ highly integrated,economical off-the-shelf data processing electronics and semiconductorlasers (laser diodes). Such features contrast with other imagingtechnologies commonly used in clinical diagnosis that require large,fixed facilities such as MRI and x-ray CT imaging. Additionally, since asignificant computational effort may be required for both imagereconstruction and data analysis, the technology particularly gains fromthe exponential growth in the ratio of computing power to cost.

It is well appreciated that optical tomography has the potential toprovide insights into anatomy and physiology that are unavailable fromother imaging methods. For example, optical tomography, using nearinfrared energy, can identify the spatial variations in blood volume andblood oxygenation levels because of its sensitivity to hemoglobinstates. These measures have considerable potential value in diagnosing abroad range of disease processes that are known to influence hemoglobinstates.

For example, a common feature of breast tumors, and solid tumorsgenerally, is the occurrence of neovascularized tissue.Ultrastructurally, these tissues are highly disorganized and exhibitfunctional abnormalities. Often the microvessels are dilated, tortuous,elongated and saccular. There is excessive branching of the vessels,including significant arterio-venous shunting as well as blind vascularendings. Aberrant vascular morphology and decreased vessel density areresponsible for increase resistance to flow. The resistance to flowcombined with an enlarged diffusion distance, due to the expansion ofthe extravascular space, can lead to perfusion with hypoxemic andnutrient-deprived blood. The net effect of this state is the occurrenceof substantial spatial and temporal heterogeneity in the tumor metabolicmicroenvironment.

Although these attributes of disease tissue are well appreciated, theavailability of a suitable detection methodology able to take fulladvantage of these characteristics has been notably lacking. Anappropriate methodology would be one sensitive not only to alteredhemoglobin states (i.e., localized variations in tissue blood volume andoxygenation states), but also to their dynamics under homeostaticconditions or in response to specific provocations.

A variety of methods involving imaging and non-imaging modalities areavailable for assessing specific features of the vasculature. Detailedimages of the vascular architecture involving larger vessels (>1 mmdia.) can be provided using x-ray enhanced contrast imaging or MRangiography. These methods however are insensitive to hemoglobin statesand only indirectly provide measures of altered blood flow. The latteris well accomplished, in the case of larger vessels, using Dopplerultrasound, and for near-surface microvessels by laser Dopplermeasurements, but each is insensitive to variations in tissue bloodvolume or blood oxygenation. Ultrasound measurements are also limited bytheir inability to penetrate bone.

In principle, imaging methods based on the detection and analysis ofscattered energy, such as optical tomograpic methods, can provide eitherdirect or indirect measures of all of these parameters. However, theknown methods and systems have several shortcomings. First, knownmethods and systems provide images having low contrast and resolution.Second, these methods and systems do not image the dynamic properties ofhighly scattering media. Third, these methods and systems requireaccurate calibration and are susceptible to errors. There are severalreasons for these problems with known systems and methods. These reasonsrelate principally to how measurements are performed and how measurementdata is analyzed.

For example, when imaging human tissue, the natural occurrence ofvascular frequencies arising from cardiac, respiratory and vasomotoractivity, produces time variations in, for example, the absorptionproperties of tissue due to changes in tissue blood volume.Significantly, the process of vasomotion, perfusion first in one region,then another, can be expected to produce spatially convolved imagesshould data be collected on a time frame that is long compared to thereciprocal of the frequency of these processes. Thus methods thatcollect time-averaged data will predictably yield images whose contrastand resolution are degraded by such variability.

Also influencing the quality of reconstructed images, is the approachused to analyze the acquired data. Many of the known imaging schemesconsider, in some manner, the comparison of measured values to predictedvalues. Typically, these methods, including that described in theBarbour '355 patent, employ numerical methods that seek to minimize thedifference between sets of measured and predicted values, and in doingso seek to provide improved estimates of the properties of the unknowntarget medium. These analysis schemes, referred to as model basedmethods, assume equivalency in the efficiency of measured detectorvalues and computed predicted values. Although the derivation ofaccurate estimates of the efficiency of measured responses is possiblein principle, in practice, the natural plasticity of tissue, its mainlyarbitrary shape and variable composition and noted variability inhemoglobin levels, all serve to confound efforts to devise practicalmethods that provide reliable estimates.

In addition, the physics of light transport in highly scattering media,such as tissue, imposes further practical constraints that relate to themethod adopted for data analysis. At issue is the required accuracy ofassumptions made in order to generate predicted detector values,especially those adopted for the initial estimate. These assumptions arecommonly referred to as the “initial guess”. Small errors in the initialguess of optical properties of the reference medium can lead to largeerrors in the computed detector values. One consequence of this can bethe severe corruption of the information content of the data vectorsleading to artifact-laden images. Adding to the mentioned uncertaintiesis the well-known property of reconstruction methods that employ linearoperators regarding their sensitivity to undetermined data sets (i.e.,insufficient amount of collected data) and measurements based onrestricted views (e.g., backscatter only, transmission only). The neteffect of these limitations can render solutions to image recoveryproblems of this type overly sensitive to the influence of experimentalnoise (i.e., ill-conditioned), provide nonunique solutions (ill-posed)or both.

These concerns are well appreciated by those skilled in the art of imagereconstruction methods. It is also well appreciated that, in general,there are no simple or well-defined methods that can be universallyapplied to overcome the noted limitations. In this regard, specificationof suitable conditions that can satisfactorily deal with the notedconcerns is an art whose successful implementation requires considerableskill.

In addition to the need to provide stable solutions to the imagereconstruction problem, consideration of the information content of thereconstructed image has considerable importance. As presently practiced,the method of optical tomography considers the evaluation of staticstates or employs time averaging methods to minimize the influence ofsignal instability originating from tissue dynamics.

The goal of these studies is to provide image maps that define spatialvariations in the optical properties of tissue (usually, absorption andscatter) from which may be derived, for instance, estimates of spatialvarying hemoglobin states. It is understood however, that the latter isfundamentally governed by dynamic processes whose details have thepotential to reveal a wealth of information regarding functionalfeatures of the vasculature, in particular as it relates to itsinteraction with surrounding tissue.

The ability to measure dynamic processes of a medium can revealinformation that is unobservable from static or time-averaged measures.In the case of physiological systems, the form of the dynamic processhas added significance. For instance, it is well understood that manytime varying processing have an underlying nonlinear character.Nonlinear dynamic processes, in biological systems, are often chaoticand exhibit the characteristic feature of sensitivity to initialconditions.

The existence of such behavior has important implications in theunderstanding of disease processes and well as for the approaches takenfor therapy. For instance, the approach needed to control a chaoticsystem is quite different from that for a linear system, wherein thesystem response is proportional to the magnitude of the input stimulus.Thus it has been proposed that more effective therapies can be realizedfrom a series of well-timed perturbations rather than from the standardapproach of applying a constant stimulus, the method commonly used inmany pharmacological interventions. Also of interest, and related tothis, is the seemingly general finding that the occurrence of chaoticbehavior in physiological systems is a sign of health and its absence isa sign of disease.

For instance, it is known that heart rate variability is chaotic.Significantly, loss of this signature with the appearance of periodicoscillations is among the strongest predictors of sudden cardiac death.A similar phenomenology has been observed in infants who succumb tosudden infant death syndrome. In this case, the normally chaoticrespiratory rate becomes periodic prior to the fatal incident.Similarly, during epileptic seizures, electroencephalographic recordingsexhibit a transition from chaotic to periodic activity.

Presently, the capacity to monitor dynamic behavior in vascularstructures is limited principally to near surface measures using laserDoppler methods. Measures of the time variability of the vascularcaliber and flow motion for larger vessels is possible using Duplexultrasound. However, these measures are insensitive to the activity ofthe microvasculature, do not provide for full cross-sectional views, andare not sensitive to the dynamics of hemoglobin states.

Although optical methods, such as Laser Doppler, pulse oximetry,photoplethysmography and the like, can be used to monitor dynamic statesof the vasulature, none are capable of providing such measures in theform of a cross-sectional image, especially in the case of large tissuestructures. Moreover, known optical methods for cross-sectional imagingof the properties of a scattering media have not been used to derivedynamic measures of these states, and are plagued with a host oftechnical limitations, such as low contrast, low resolution images,prolonged computing times, excessive sensitivity to errors in initialestimates. These limitations render it unlikely that useful informationregarding dynamic states could be derived from these known methods.

Overcoming the indicated drawbacks and concerns is critical forwidespread practical implementation of optical tomography as adiagnostic tool because (1) improved contrast and resolution areessential to feature identification and visualization, (2) a static(snapshot) or time averaged image of a time evolving property does notprovide discovery of the physiological dynamic processes and (3)measures of dynamic processes can yield critical information needed forimproved diagnostic methods and therapies.

For the foregoing reasons, there is a need for a method of improving thecontrast and resolution of reconstructed images. There is also a needfor a method of imaging dynamic properties of dense scattering media,especially as it relates to dynamic properties of vascular states inlarge tissue structures as revealed by time variation in hemoglobinstates. There is yet a further need for a method that can providedynamic images without undue reliance on complex calibration schemes orcomputationally intensive numerical methods.

SUMMARY OF THE INVENTION

The present invention satisfies these needs by providing a method for(1) enhancing resolution and contrast of various dynamic features of themedium being imaged, (2) imaging the dynamic properties of the targetmedium, and (3) generating images of the dynamic properties usingtechniques that reduce the need for undue system calibration and producemore stable solutions with the reconstruction algorithm.

It is one object of the invention to generate a map of the dynamicproperties of a scattering medium in a cross-sectional view from a timeseries of collected data measurements. The measurements are obtained bydirecting energy into the target during a period of time and measuringthe energy emerging from the target during the period of time, whereby atime series of measured energy is collected.

It is another aspect of the invention to generate a map of thecross-sectional dynamic properties of the target medium in across-sectional view from a time series of images of the cross-sectionalproperties of the target medium, wherein the dynamic properties areextracted from the time series of images using time series analysismethods.

It is a further aspect of the invention to generate a map of thecross-sectional dynamic properties of the target medium from the timeseries of measurements at each detector, wherein the dynamic propertiesare extracted from the time series of measurements at each detectorusing time series analysis methods.

It is a further aspect of the invention to construct a map of thedynamic properties of the medium, wherein the measured energy isprocessed using a modified perturbation formulation referred to as thenormalized difference method, based on the radiation transport equation,using relative energy measurements. The relative energy measurementsbeing, for example, the relative difference between an instantaneousmeasure and a time average mean of measures.

BRIEF DESCRIPTION OF THE FIGURES

For a better understanding of the invention, together with the variousfeatures and advantages thereof, reference should be made to thefollowing detailed description of the preferred embodiments and to theaccompanying drawings wherein:

FIG. 1 is a schematic illustration of an exemplary imaging system;

FIG. 2A is an illustration of a mechanism used to rhythmically inflatetwo balloons;

FIG. 2B is a plan view of the illustration of FIG. 2A;

FIG. 2C is an illustration of the geometric arrangement of opticalfibers;

FIG. 3 is a reconstructed image of the cross-sectional properties of atarget medium shown in FIG. 2A at an instant in time;

FIG. 4 is a cross-sectional map of the 0.2 Hz component of the amplitudeof the image sequence's DFT for the target medium shown in FIG. 2;

FIG. 5 is a cross-sectional map of the 0.2 Hz component of the phase ofthe image sequence's DFT for the target medium shown in FIG. 2 with nophase difference between the oscillating balloons;

FIG. 6 is a cross-sectional map of the 0.2 Hz component of the phase ofthe image sequence's DFT, 180° phase difference between oscillatingballoons;

FIG. 7A is a schematic illustration of apparatus with two-balloons;

FIG. 7B is a series of reconstructed cross-sectional images from DFTanalysis of image time series at 0.1 Hz and 0.24 Hz;

FIG. 8 is a cross-sectional map of ratio of FT amplitude of cardiac torespiratory frequency from forearm study;

FIG. 9 is an reconstructed MR image of forearm, with identifiedanatomical structures;

FIG. 10 is an overlay of images shown in FIGS. 8 and 9;

FIG. 11A is a graph of the CSD spectrum amplitude at position (12, 19);

FIG. 11B is a graph of the CSD spectrum amplitude at position (11, 10);

FIG. 11C is a graph of the CSD spectrum amplitude at position (18, 27);

FIG. 12 is a cross-sectional map of amplitude, ×10⁴, of the DFT at thefinger-flex frequency;

FIG. 13 is an overlay of the map shown in FIG. 12 and an anatomicalimage (MRI) of the same tissue structures;

FIG. 14 is a graph of the temporal variations in μ_(a) (absorptioncoefficient) for pixels located within the involved muscle groups;

FIG. 15 is a graph of the time-dependent intensity variations forforearm measurements observed for a detector positions opposite thesource;

FIG. 16 is a cross-sectional map of the correlation dimension, v,computed, in each pixel, from time series of reconstructed images ofμ_(a) (810 nm) obtained from representative detector data shown in FIG.15;

FIG. 17 is a cross-sectional map of the correlation dimension, v,computed in each pixel from time series of reconstructed images of μ_(a)(810 nm) obtained from measurements performed on a solid white Delrin®rod;

FIG. 18 is a cross-sectional map of λ₁ (largest value of Lyapunovexponent) in each pixel derived from the same time series ofreconstructed μ_(a) (810 nm) images as used to generate v map shown inFIG. 16; and

FIG. 19 is a cross-sectional comparison map of the λ₁ map in FIG. 18 tocorresponding results obtained from surrogate data sets.

DETAILED DESCRIPTION OF THE INVENTION

Introduction

The present invention is built upon known basic approaches totomographic imaging of a scattering medium, by enabling theinvestigation of the spatio temporal dynamics of a target medium. Thisenhancement represents a fundamentally new imaging modality with broadbased applications such as in industrial processes involving preparationof powdered mixtures, imaging subsurface turbulence in murky waters andthroughout clinical medicine both as a highly sensitive diagnosticimaging tool, as well as for use in treatment planning, therapeuticmonitoring and follow-up.

For example, the present invention recognizes that the occurrence ofdynamic behavior in human tissue lends the collected data andreconstructed images to evaluation using time series analysis methods,including pattern recognition schemes. Dynamic behavior is known tooccur in vascular structures in response to cardiac, respiratory,vasomotor and local metabolic influences. Because optical tomography,based on near infrared measurements is sensitive to hemoglobin, which isnearly always restricted to the vascular compartment time variations inthe measured signal, can be used to reveal the vascular frequencies.Although, similar information can be obtained by photoplethysmography,it is not available for large tissue structures nor in ancross-sectional imaging modality.

Detection of these temporal signatures/dynamic behaviors in across-sectional view provides a wealth of new information. For instance,the presence (or absence) of a cardiac frequency with respect to knownanatomic landmarks (e.g., major arteries) provides a measure of thepatency of localized arterial flow. Low-amplitude signals might suggestthe presence of arterial stenosis. In addition, it is well appreciatedthat blood flow is under autonomic control. The absence or attenuationin amplitude of the respiratory or vasomotor frequencies might suggestan underlying peripheral neuropathy.

Further insight concerning the presence or absence of basic controlmechanisms can be obtained from time-frequency analysis methods. Forinstance, the low-frequency vasomotor response in the major arteries ofthe forearm are nearly synchronous, both ipsilaterally andcontralaterally, and under tight autonomic control. Disorganization ofthis control might suggest the presence of pathology. In fact, loss ofcontrol mechanisms is a well known hallmark of neoplastic tissues. Thenearly total absence of the vasomotor response in tumor tissue is ofparticular interest.

Still another factor motivating dynamic measures is an appreciation ofthe information content available from the study of homeostaticchallenges. It is well appreciated that the altered vasculararchitecture in tumors leads to a condition of sluggish perfusion. Theexistence of such states maybe be clearly revealed by induction of ahomeostatic challenge. This could take on several forms. For example,simple manipulations such as deep breathing, response to cold, or mildvenous congestion imposed by inflation of a pressure cuff, all have beenshown to produce large amplitude local variations in the optical signalsattributable to tissue blood volume in the forearm. Extension of thesemeasures to the breast and other tissue could reveal the presence oftumors with high specificity by virtue of an altered hemodynamicresponse.

There are other important implications to be considered concerning thevalue of measuring the time-varying properties of hemoglobin states in across-sectional view. For instance, the details of the measured dynamicsat the various vascular frequencies need not be the same for measures ofblood volume and blood oxygenation. In addition, given the knownheterogeneity in tissue perfusion, even in healthy tissue, whateverdynamic features do exist will themselves be spatially varying.Moreover, altered dynamics in blood volume or blood oxygenation indicatealtered local metabolic states, the existence of more central controldeficits (e.g., autonomic, cardiac, respiratory), or both. Hence thesemeasures can provide an assessment of an integrated physiological state,and also possess important features pertaining to tissue/vascularcoupling. In fact, all of the spatiotemporal features attributable tohemoglobin states that do exist will likely respond to a host ofpharmacological agents and other treatment modalities. Such measurescould serve to identify desired or undesired responses to therapy.

Although these examples relate to the measure of hemoglobin states, itwill be appreciated by one skilled in the art that there are a plethoraof dynamic features that may be observable using tomographic measures inscattering media in human tissue as well as other scattering media, suchas murky water, foggy atmospheres and the like. For instance nerveactivation and muscle movement can also by measured owing to changes inscattering properties of tissue. In addition, various exogenous agentsthat produce contrast in absorption, scatter or fluorescence can beintroduced to observe dynamic states of tissue.

The present invention includes new data collection schemes and imageanalysis methods from which are extracted maps that reveal spatialtemporal signatures of the target medium. These methods can include,linear and non-linear time series analysis methods as well as techniquesfor pattern recognition.

There are three principal elements to practical dynamic imaging. Thefirst element is the use of a fast, parallel, multi-channel acquisitionsystem. For example, a system for use in dynamic optical tomographicimaging of tissue is described briefly below and in further detail inthe copending Barbour 4147PC1 application. The second element is toevaluate the acquired tomographic data to produce a cross-sectionalimage that reveals perturbations of the reconstructed opticalcoefficients from, for example, a time-averaged mean. This perturbationmethod is described briefly below and in further detail in the Barbour4149PC2 application. The third element is to collect a time series ofdata and subject either the time series of measured data or a timeseries of reconstructed images from the measure data to analysis. Thiscan include use of various linear and nonlinear time-series analysismethods, and related techniques such as pattern recognition, to extractfeatures identifying dynamic behavior in the properties of the targetmedium. These methods are discussed in further detail below.

The System

Numerous imaging systems such as those disclosed in the Barbour '355patent, and the Barbour '322 patent have been developed for use inimaging of a scattering medium which are hereby incorporated byreference.

A system providing high speed data capture of one or more wavelengthssimultaneous is disclosed in the Barbour 4147PC1 application. Thissystem is capable of capturing multiple wavelength data at rates up to150 Hz and enables the reconstruction of cross-sectional images ofreal-time events associated with vascular reactivity in a variety oftissue structures (e.g., limbs, breast, head and neck). Fast datacollection methods are particularly useful because there are manydisease states with specific influences on the spatial-dynamicproperties of vascular responses in hemoglobin states.

A schematic illustration of an exemplary optical system for fast datacollection is shown in FIG. 1. This system includes a computer 102,energy sources 104, 106, a source demultiplexer 108, an imaging head110, detectors 112 and a data acquisition board 114.

A target 116 placed in the imaging head 110 is exposed to optical energyfrom sources 104, 106. The optical energy originating from energysources 104, 106, is combined by beam splitter 118 and is delivered tosource demultiplexer 108. The source demultiplexer 108 is controlled bycomputer 102 to direct the optical energy to source fibers 120sequentially. Although two energy sources 104, 106 are shown in thisembodiment, an additional number of energy sources, each havingdifferent wavelengths can be employed. Moreover a single variablewavelength energy source can be implemented such as Ti-Sapphire laser,dye lasers and the like.

Each source fiber 120 carries the optical energy from the demultiplexer108 to the imaging head 110 where the optical energy is directed intothe target 116. The imaging head 110 contains a plurality of sourcefibers 120 and detector fibers 122 for transmitting and receiving lightenergy, respectively. Each source fiber 120 forms a source detector pairwith each detector fiber 122 in the imaging head 110 to create aplurality of source detector pairs. The optical energy entering thetarget 116 at one location is scattered and may emerge at any locationaround the target 116. The emerging optical energy is collected bydetector fibers 122 mounted in the imaging head 110.

The detector fibers 122 carry the emerging energy to detectors 112. Thedetectors 112 measure the energy density of the collected optical energyand generate a corresponding signal. The detectors may comprise anyknown photodetector such as avalanche photodiodes (APD), silicon PINphotodiodes, silicon photodiodes (SiPD), charged couple devices (CCD),charge inductive device (CID), photomulitplier tubes (PMT),multi-channel plate (MCP) and the like. The energy density is, forsituations in which propagating energy is any type of electromagneticradiation, equivalent to the intensity of said radiation.

The data acquisition board 114 receives the signal, separates it bywavelength, and samples and holds the separated signals for deliver tocomputer 102. The computer 102 in turn reads and stores the signal foruse in image reconstruction and other analysis.

Although the above described systems uses DC measurement techniques withan near infrared energy source, it will be appreciated by those skilledin the art that similar systems can be implemented using othermeasurement techniques, such as time resolved measurements and frequencydomain methods. In addition measures of acoustic signals produced inresponse to optical absorption (i.e., the photo acoustic effect) as wellis use of acoustic sources can be similarly implemented.

Data Collection, Data Analysis and Image Reconstruction In prior methodsfor imaging of scattering media, measured data was collected for atarget medium and an image of the cross-sectional properties of thetarget was generated based on the measured data. The present inventionimproves on this basic process by (1) collecting a time series ofmeasured data and (2) analyzing either the raw measured data or a timeseries of reconstructed images to extract information relating to thedynamics of the target medium.Data Collection

As discussed above, previously known systems and techniques sought tocollect either a time averaged measure thereby minimizing the effect ofdynamic behavior or to collect a snapshot of the target. Both of thesetechniques fail to recognize the value of information contained in thetime domain.

The present invention recognizes the value in measuring dynamicproperties in highly scattering media and provides methods forcollecting and processing the measured data to generate images ofdynamic properties. These methods rely on the high speed collection ofdata using a system such as that described above so that a time seriesof data sets can be collected in the time domain. The time seriesanalysis techniques, discussed further below, can then be used toextract the dynamic properties of the target from either the measuredraw data directly or from a time series of images reconstructed from themeasured data.

Data is collected in the present invention using high speed imagingsystems such as the system illustrated in FIG. 1. By way of example,referring to FIG. 1, a measured data set is obtained by deliveringoptical energy to each of the source fibers 120 sequentially andmeasuring the emerging optical energy at each detector fiber 122 inparallel for each sequential source location. This process is repeatedover a period of time so that a time series of data sets are generated,each data set represent a complete set of data (e.g., detector readingfor each source position) for reconstruction of a cross-sectional imageat an instant in time.

Reconstruction

As discussed above, the time series analysis techniques may be employedprior to image reconstruction (i.e., on the raw measured data) or afterimage reconstruction (i.e., on the time series of reconstructed images).Where time series analysis methods are applied to raw measured data,image reconstruction will follow based on the processed and extracteddynamic information. For example, one method may compute the discreteFourier transform of the measured data and then reconstruct images basedon the Fourier coefficients at, for instance, a selected frequency.

Of the known reconstruction methods, one group of methods often adoptedto analyze the measured data are perturbation methods. Exemplarytechniques are disclosed in detail in the Barbour '355 patent and theBarbour 4149PC2 application. Briefly, in model-based techniques aperturbation formulation based on the radiation transport equation, orits approximation the diffusion equation, is generated. The radiationtransport equation is a mathematical expression describing thepropagation of a particle through a medium from a source location to adetector location as a function of the properties of the medium. Theperturbation formulation of this equation relates a perturbation in themeasure energy density at a detector to a corresponding perturbation inthe properties of the medium.

The standard perturbation equation has the following forms:W·δx=δ1  (1)

In equation (1), δ1 is a vector of source-detector pair intensitydifferences between the measured target medium and the known referencemedium (i.e., δI=I−I_(r)). W is the weight matrix describing theinfluence that each volume element (“voxel”) of the reference medium hason energy traveling from each source to each detector, for allsource-detector pairs. The volume elements are formed by dividing aslice of the reference medium into an imaginary grid of contiguous,non-overlapping pieces. Physically, the weight matrix contains the firstorder partial derivatives of the detector responses with respect to theoptical coefficients of each volume element of the reference medium. δxis the vector of differences between the known optical properties (e.g.,absorption and scattering coefficients) of each volume element of thereference medium and the corresponding unknown optical properties ofeach volume element of the target medium.

A modification to this standard perturbation equation is presented anddisclosed in the Barbour 4149PC2 application. This modified perturbationequation employs relative measurements of energy at the target detectors(e.g., the relative difference between an instantaneous measure and atime averaged mean of measures at a detector). This method improves thestability of the solution to the perturbation equation and reducessensitivity of the equation to the selected reference medium.

The modified standard perturbation equation replaces δI in the standardperturbation equation with a proportionate relative difference betweentwo measured values multiplied by a reference term of the required unitsas set forth in the equation (2) below: $\begin{matrix}{{\delta\; I_{r}} = {\left( \frac{I - I_{0}}{I_{0}} \right) \cdot I_{r}}} & (2)\end{matrix}$

In equation (2), I_(r) is the computed detector reading corresponding toa source-detector pair of a selected reference medium, and I and I₀represent two data measurements for a corresponding source-detector pairon one or more targets (e.g., background vs. target, or time-averagedmean vs. a specific time point, etc.). The resultant term δI_(r)therefore represents a proportionate relative difference between twosets of measured data This modified equation limits the effects ofmodeling errors and minimize ill-conditioning of the inverse problemwhile retaining the correct units in the solution.

The corresponding perturbation equation using this modified term is:W _(r)˜δx=δ1_(r)  (3)

In equation (3) W_(r) and δx are the same as W and δx in equation (1).Referring to equations (2) and (3), it can be seen that in the limit,when I_(r) equals to I₀, this equation reduces to the standardperturbation formulation shown in equation (1). This formulationrepresents the Born approximation formulation of the modifiedperturbation equation. A similar expression may be written for the Rytovapproximation in the following form: $\begin{matrix}{{l_{n}\left( \frac{I}{I_{0}} \right)} = {\frac{W_{r}}{I_{r}}\delta\; X}} & (4)\end{matrix}$

The inventive operation performed by equation (2) amounts to aprojection procedure by which the information content of a measuredproportionate relative data vector is projected onto a selectedreference medium. Apart from simplifying measurement requirements, themethod represented by equations (3) and (4) also reduces susceptibilityof the perturbation equation to boundary and optical property variationbetween the target and the reference medium. Additionally, iterativesolutions of equations (3) and (4) can be easily implemented. Inprinciple, the techniques used in the modified perturbation equation canbe used to evaluate any measured proportionate relative quantity to areference quantity, for example acoustic signals.

Time Series Analysis

As discussed above, time series analysis methods may be applied toeither the time series of measurement data or time series ofreconstructed images. Applying these methods, the present inventionenables investigation of the temporal spatial dynamics of a targetmedium.

Time series analysis techniques are well-known to provide measures ofthe frequency structure and time correlation of time varying processes.In addition, hybrid methods such as time frequency analysis (e.g., shorttime Fourier transform, wavelet analysis) can be used to reveal morecomplex behaviors in physiological systems such as the known frequencymodulation of the cardiac rate upon respiration. Other linear methodscan be applied to time domain data such as principle component analysisto reveal the time evolution of specific processes in image maps.Practical examples of these methods as shown in FIGS. 2–19 can be usedto identify, in the case of frequency analysis, specific features of thevascular tree, in particular, in the case of imaging human tissue, themajor arteries in the forearm. This is possible because of the knownstructure dependent frequency response of the vasculature (e.g., majorarteries beat at a cardiac frequency). Further, time correlation methodsmay be employed to isolate anatomical landmarks having known phaserelationships (e.g., action of antagonistic muscle groups).Additionally, similar procedures can be used on the resultant spatialmap to produce the relevant orthogonal information associated with thesefunctions (i.e., maps of the amplitude or phase vs. frequency,correlation vs. time lag, etc.). In addition, complex behaviors in theform of nonlinear chaotic activity is often observable in measures ofphysiological systems. As is shown in FIGS. 15–19, the describedinvention is capable of detecting such signatures in a cross-sectionalview, in particular as it relates to chaotic behavior of the vasculatureas measured by time variations in hemoglobin states. While such behaviorhas been observed in surface measurements of tissue, with the exceptionof the instant invention, no other imaging technique is capable ofgenerating equivalent information.

Experimental Validation

The following discussion and some of the above examples presents resultsvalidating the methods and advantages of the present inventions. Theseexamples are presented merely as an illustration of the benefits of thepresent invention and in particular of the ability of dynamic opticaltomography to provide high contrast images of time varying features indense scattering media.

Contrast and Resolution Enhancement

As disclosed above, application of these time-series analysis methodscan extract significantly greater contrast and resolution from theoptical imaging data. Image resolution and contrast are importantfactors for any imaging method. Experience with optical tomography hasshown that recovered images have low resolution and contrast. As anexample, referring to FIG. 2, tomographic data was collected from alaboratory vessel measuring approximately 8 cm in diameter andcontaining four balloons 10-1 through 10-3, two of which are static andthe other two were modulated at a fixed frequency, with a relative phaseof either 0° or 180°. The balloons 10 were suspended in a backgroundscattering medium 12 consisting of 2% (v/v) Intralipid and were filledwith a dilute hemoglobin solution (50 μM). The balloons 10 were arrangedin the vessel 14 as shown in FIGS. 2A and 2B. An example of areconstructed cross-sectional image of absorption obtained at a singletime point is shown in FIG. 3. Inspection of FIG. 3 shows that whereasthe four objects structure is detectable, the individual balloons arenot well resolved, and the recovered contrast is poor. In contrast tothis image quality, the map shown in FIG. 4, revealing the discreteFourier transform (DFT) amplitude at the 0.2 Hz modulation frequency hasnearly 100 fold increased contrast. Also evident is that the spatialresolution of the DFT map is considerably improved compared to the mapsshown in FIG. 3. The phase portions of the DFT's at 0.2 Hz for thein-phase and the 180°-out-of-phase cases are shown in FIGS. 5 and 6. Thephase relations between the two modulated balloons is correctlyrecovered in each case.

A similar experiment resolving internal structures based on differencesin modulation frequency is shown in FIGS. 7A and 7B. In FIG. 7A twoballoons were present, one modulated at 0.1 Hz, the other at 0.24 Hz.Inspection of the reconstructed images shown in FIG. 7B, at thecorresponding modulation frequencies shows complete spatio-temporalresolution.

Dynamic Physiological Imaging Studies on the Human Forearm

Following a protocol similar to that used to investigate dynamicbehavior in the laboratory vessel, the information retrievable fromdynamic studies performed on the human forearm at rest and in responseto specific provocations are also explored.

Arm at Rest

The natural occurrence of vascular frequencies due to respiratory andcardiac activity can be exploited to produce a spatial map revealing thepresence of different functional components of the vascular tree. FIG. 8shows a map of the logarithm of the ratio of the computed DFT amplitudesat the cardiac and respiratory frequencies obtained from a times-seriesof measurements on the forearm at rest. FIG. 9 is a representative MRimage in the same region of the forearm. An overlay of the two mapshaving the same orientation is shown in FIG. 10. Inspection reveals thatin the vicinity of the radial 1, interosseous artery 3 and ulnar artery5, the ratio of the DFT amplitudes (cardiac to respiratory) is nearlyten times larger than it is in other regions. The radius 2 and ulna 4are also shown. This response can be seen more clearly in FIGS. 11A–11C,which shows the amplitude of the cross-spectral density (CSD) between asurface detector and specific locations in the image. The particularspectra shown were obtained from points in the image corresponding tolocations in the flexor digitorum superficialis muscle, and points nearthe radial and interosseous arteries. Inspection clearly reveals that inmuscle, which contains a dense matrix of microvascular structures, thedominant signal coincides with the respiratory frequency. In contrast,we observe a dominant cardiac frequency in the vicinity of the majorarteries. These findings coincide well with the known hemodynamicbehavior of the different structures in the vascular tree.

Finger-Flex Study

In this study we explored further our ability to measure dynamicbehavior in highly scattering media, such as the human forearm, byexamining a time series of images obtained from measurements performedon a subject while conducting a finger-flex exercise. Finger flexinginvolves the action of so-called antagonistic muscle groups that arelocated on opposite sides of the forearm, specifically, the flexordigitorium superficialis on the ventral side and the extensor digitorumon the dorsal side. Results shown in FIG. 12 show a map of the amplitudeof the DFT at the finger-flexing frequency (˜0.25 Hz). FIG. 13 shows anoverlay of this image onto an MR image of the same forearm in the sameorientation. Inspection reveals that positions of maximum amplitude forfinger-flexing coincide well with the two involved muscle groups.Further evidence supporting the accuracy of this assignment is shown inFIG. 14. Plotted are time-series values of recovered absorption valuesat points in the image coinciding with the involved muscles. Noteworthyis the observation that the two signals are approximately 180° out ofphase with each other, which is the expected response from the action ofantagonistic muscle groups.

Imaging the Complexity of the Vascular Response

One feature of the disclosed methods that may have considerablepractical significance is the potential to obtain reliable estimates ofthe complexity of the vascular response. The significance of suchmeasures is that they can reveal features of an integrated physiologicalstate that are simply not measurable by using linear methods.Accordingly, in this study we have sought to confirm previous reports inthe literature that indicate nonlinear chaotic behavior in vascularrhythms, as has been measured typically by surface laser Dopplermethods. Specifically, we have computed certain measures that can revealnonlinear behavior from time-series image data of the forearm of asubject performing a series of deep-breathing exercises. The influenceof a respiratory stimulus was chosen as a simple noninvasive means toamplify the natural oscillatory activity of the vasculature. FIG. 15shows an example of the magnitude of the signal response caused by deepbreathing. The measures we have computed from the pixel data include v,the correlation dimension, and, λ₁, the largest Lyapunov exponent.Positive values for the latter indicate divergence of trajectories instate space, a finding consistent with a chaotic time series. A value ofzero indicates no net divergence and is observed with stochastic data,while negative values indicate system damping.

The computed values displayed in a cross-sectional shown in FIG. 16,below were obtained using the method of Rosenstein et al. (Physica D,65, 117–134 (1993)) for small data sets. The image series comprised 240consecutive time points reconstructed from data collected at a rate of2.8 Hz. The respiration rate was 0.08 Hz for the first 150 time points,and then doubled for the remainder of the measurement period. To assistin timing and repeatability, the volunteer was asked to adjust hisrespiratory rate to match the beat of a metronome.

Referring to FIG. 16, each pixel in the derived image series wasband-pass filtered over a range that included the vasomotor andrespiratory frequencies (0.05–0.35 Hz). Such filtering colors the noiseand can lead to spurious low values for the correlation dimension. Forthis reason data obtained from a stochastic time series in a similarfashion was treated similarly. The control data set was obtained bycomputing an image time series from measurements performed on a solidwhite plastic rod composed of white Delrin®, whose scattering propertiesare similar to those of tissue. The rod diameter was 9 cm.

As shown in FIG. 16, the range of correlation dimensions computed fromthe physiological image time series is between 2 and 4 for most pixels.Significantly, these values are in close agreement with previous studiesthat measured the correlation dimensions of arterial and arteriolarvasomotion in in vitro rat or rabbit preparations. A similar analysis onthe rod data, shown in FIG. 17, reveals significantly higher values, afinding consistent with the character of the time series.

The image in FIG. 18 shows a spatial map of λ₁ computed for each pixel.Because a stochastic time series can also yield positive values for λ₁,we have computed surrogate data sets for each pixel and performed a testfor statistical significance. Surrogate data were computed byrandomizing the phase of the Fourier spectrum followed by recovery ofthe resultant time series from which the surrogate λ₁ were computed.FIG. 19 shows a map of the significance level, in each pixel, for thephysiological data compared to surrogate data. Inspection reveals thatmost areas of the map have positive λ₁ values that fail the nullhypothesis. Note that the optical measures performed are known toexhibit considerable sensitivity to the microvessels. These findings,indicating the occurrence of nonlinear chaotic activity of thevasculature revealed in a cross-sectional view, thus coincide well withprevious reports whose measures were limited to the study of vasculardynamics in the near surface of tissue. Its worth emphasizing that theabove findings are the first of their kind and, to the best of ourknowledge, no other imaging modality is capable of providing equivalentmeasures of vascular dynamics.

Strategies of the types explored can be used to investigate a range ofvascular anomalies. For example, the deep breathing exercise may aid inthe diagnosis of cancer. Solid tumors have altered vascular beds.Enhanced breathing will cause redistribution of the flow patterns intissue. Regions having altered vascular architecture should behave verydifferent from those with a normal architecture. A similar approach canbe used to identify tissue regions “at risk” due to the effects ofdiabetes, atheroscleorsis (small vessel disease), or long term effectsof smoking (large vessel disease).

Although the examples above focus on near infrared energy sources forimaging human tissue, the methodology of the present invention isapplicable to essentially any wavelength for any energy source (e.g.,electromagnetic, acoustic, etc.), responsive to dynamic properties of ascattering medium and for any source condition (e.g., time independent,time harmonic, time resolved). Moreover, although the examples aboveinvestigate cross-sectional imaging, it will be appreciated that thesemethods are applicable to and include volumetric measures and imaging.

Further, although illustrative embodiments have been described herein indetail, those skilled in the art will appreciate that variations may bemade without departing from the spirit and scope of this invention.Moreover, unless otherwise specifically stated, the terms andexpressions used herein are terms of description and not terms oflimitation, and are not intended to exclude any equivalents of thesystem and methods set forth in the following claims.

1. A method for enhanced imaging of a target medium comprising: usingoptical tomography to direct energy into a target medium from at leastone source during a period of time, the target medium having dynamicproperties during the period of time; wherein the energy from the atleast one source is highly scattered by the target medium and emergesfrom the target medium at different locations around the target medium;measuring the energy emerging from the target medium during the periodof time using a plurality of detectors positioned to detect the emergingenergy at the different locations; the energy emerging from the targetmedium at the different locations being a function of the dynamicproperties of the target medium; and generating a map of the dynamicproperties of the target medium based on the measured energy emergingfrom the target medium at the different locations.
 2. The method ofclaim 1, further comprising: generating a time series of images of theproperties of the target medium based on the measured energy emergingfrom the target medium, wherein each image represents thecross-sectional properties of the target medium at a time intervalduring the period of time.
 3. The method of claim 2, wherein generatingthe map of the dynamic properties of the target medium comprisesprocessing the time series of images using time series analysis methods.4. The method of claim 1, wherein generating the map of the dynamicproperties of the target medium comprises processing the measured energyat each of the plurality of detectors using time series analysismethods.
 5. The method of claim 1, wherein the map of the dynamicproperties is generated using time series analysis methods that compriselinear time series methods.
 6. The method of claim 5, wherein the lineartime series analysis methods are at least one of frequency analysis,time correlation analysis, time frequency analysis and principlecomponent analysis.
 7. The method of claim 1, further comprisingapplying a provocation to the target medium.
 8. The method of claim 7,wherein the target medium comprises human tissue having a vascular treeand the provocation has a dynamic effect on the vascular tree.
 9. Themethod of claim 8, wherein the provocation comprises an autonomicstimulus.
 10. The method of claim 8, wherein the provocation comprises alocal stimulus.
 11. The method of claim 1, wherein the energy comprisesoptical energy having a wavelength in the near infrared region of theelectromagnetic spectrum.
 12. The method of claim 11, wherein theoptical energy directed toward the medium includes at least twowavelengths of near infrared energy.
 13. The method of claim 11, whereinthe target medium comprises human tissue having a vascular treecontaining blood, the vascular tree comprising veins, arteries and microvessels, the blood having time varying absorption and scatteringproperties to the near infrared energy as a function of bloodoxygenation and blood volume.
 14. The method of claim 13, whereingenerating a map of the dynamic properties of the target mediumcomprises generating an image of at least one of the time varyingabsorption and scattering properties of the target medium.
 15. Themethod of claim 13, wherein generating a map of the dynamic propertiesof the target medium further comprises using time series analysis toenhance the contrast of at least one of veins, arteries and microvessels.
 16. The method of claim 1, wherein the energy comprises opticalenergy in the visible spectrum.
 17. A system for enhanced imaging of atarget medium, comprising: a source for using optical tomography todirect energy into a target medium from at least one source during aperiod of time, the target medium having dynamic properties during theperiod of time; wherein the energy from the at least one source ishighly scattered by the target medium and emerges from the target mediumat different locations around the target medium; a plurality ofdetectors positioned at the different locations for measuring the energyemerging from the different locations of the target medium during theperiod of time, the energy emerging from the target medium at thedifferent locations being a function of the dynamic properties of thetarget medium; a data acquisition means for receiving the measuredenergy during the period of time; and a computer responsive to the dataacquisition means, the computer having code for generating a map of thedynamic properties of the target medium based on the measured energyemerging from the target medium at the different locations.
 18. Thesystem of claim 17, wherein the computer further includes code forgenerating a time series of images of the properties of the targetmedium based on the measured energy emerging from the target medium,wherein each image represents the cross-sectional properties of thetarget medium at a time interval during the period of time.
 19. Thesystem of claim 18, wherein the computer further includes code forprocessing the time series of images using time series analysis methodsto generate the map of the dynamic properties.
 20. The system of claim19, wherein the computer further includes code for processing themeasured energy at each detector using time series analysis methods togenerate the map of the dynamic properties of the medium.
 21. A methodfor enhanced imaging of a target medium comprising: directing energyinto a target medium from at least one source during a period of time,the target medium having dynamic properties during the period of time;measuring the energy emerging from the target medium during the periodof time using at least one detector, the energy emerging from the targetmedium being a function of the dynamic properties of the target medium;and generating a map of the dynamic properties of the target mediumbased on the measured energy emerging from the target medium, and usingnonlinear time series analysis methods.
 22. A method for enhancedimaging of a target medium comprising: directing energy into a targetmedium from at least one source during a period of time, the targetmedium having dynamic properties during the period of time; measuringthe energy emerging from the target medium during the period of timeusing at least one detector, the energy emerging from the target mediumbeing a function of the dynamic properties of the target medium; andgenerating a map of the dynamic properties of the target medium based onthe measured energy emerging from the target medium; wherein the energyis measured by a plurality of measurements collected at a sampling ratenot less than twice the reciprocal of the highest frequency of a dynamicproperty to be imaged.
 23. A method for enhanced imaging of a targetmedium comprising: directing energy into a target medium from at leastone source during a period of time, the target medium having dynamicproperties during the period of time; measuring the energy emerging fromthe target medium during the period of time using at least one detector,the energy emerging from the target medium being a function of thedynamic properties of the target medium; and generating a map of thedynamic properties of the target medium based on the measured energyemerging from the target medium; wherein: the measured energy isprocessed using a modified perturbation formulation of a radiationtransport equation; and the modified perturbation formulation usesrelative energy measurements.
 24. The method of claim 23, wherein therelative energy measurements are the relative differences between ameasure at an instant in time and a time average mean of measures duringthe period of time.